{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimate average velocity from displacement time-series\n",
    "\n",
    "Given the displacement time-series in meters $d^i, i=1,...,N$, with $d^i = -\\frac{\\lambda}{4\\pi}\\phi_{dis}^i$, the average velocity is estimated as:\n",
    "\n",
    "$d^i = vt_i + c, i=1,...,N$\n",
    "\n",
    "The standard deviation of the estimated velocity is given by equation (10) in Fattahi and Amelung (2015, JGR) as:\n",
    "\n",
    "$\\sigma_v = \\sqrt{\\frac{\\sum_{i=1}^N (d_i - \\hat d_i)^2}{(N-2)\\sum_{i=1}^N (t_i - \\bar t)^2}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Displacement [m]: [0.         0.01039512 0.08327228 0.05231298 0.08403516 0.02249366\n",
      " 0.05272516 0.06046131 0.04574046 0.02036186 0.04161016 0.07467044\n",
      " 0.09952251 0.08487235 0.05968709 0.05444766 0.11809858 0.07724472\n",
      " 0.12519749 0.06858165 0.1234072  0.11500141 0.09980721 0.12312403\n",
      " 0.07900213 0.15446893 0.08542804 0.16449594 0.10920841 0.08936961\n",
      " 0.14076454 0.08593687 0.16099531 0.15815976 0.15378566 0.09737901\n",
      " 0.13644764 0.12331895 0.16866511 0.19502585 0.15885452 0.11710548\n",
      " 0.13044944 0.18019792 0.19854746 0.13466152 0.17910256 0.22488263\n",
      " 0.19199365 0.23179678]\n",
      "Dates [year]: [2014.5        2014.53285421 2014.56570842 2014.59856263 2014.63141684\n",
      " 2014.66427105 2014.69712526 2014.72997947 2014.76283368 2014.79568789\n",
      " 2014.82854209 2014.8613963  2014.89425051 2014.92710472 2014.95995893\n",
      " 2014.99281314 2015.02566735 2015.05852156 2015.09137577 2015.12422998\n",
      " 2015.15708419 2015.1899384  2015.22279261 2015.25564682 2015.28850103\n",
      " 2015.32135524 2015.35420945 2015.38706366 2015.41991786 2015.45277207\n",
      " 2015.48562628 2015.51848049 2015.5513347  2015.58418891 2015.61704312\n",
      " 2015.64989733 2015.68275154 2015.71560575 2015.74845996 2015.78131417\n",
      " 2015.81416838 2015.84702259 2015.8798768  2015.91273101 2015.94558522\n",
      " 2015.97843943 2016.01129363 2016.04414784 2016.07700205 2016.10985626]\n",
      "True velocity [m/year]: 0.1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAUV0lEQVR4nO3df6wlZXnA8e+zu6xW1IgLKgLLQrrWbFto4Qa2UeuvSME2XVv/QahaW7IlkVhiaYv9lVpNW432V6SuG0qjCYbEKuk2pYIxWqOy7d5VRBZcXFdWFmhZFuLPFFj36R9nbnu9nHvvnHvnnJl5z/eTnOz5MTPnPefOPuedZ55538hMJEnlWtN2AyRJ42Wgl6TCGeglqXAGekkqnIFekgq3ru0GDHPyySfnpk2b2m6GJPXG3r17H8nMU4a91slAv2nTJmZnZ9tuhiT1RkQcWuw1UzeSVDgDvSQVzkAvSYUz0EtS4Qz0klQ4A70kFc5AL0kdsPfQY1z3mQPsPfRY49vuZB29JE2TvYce4/Lrd/PEseOsX7eGG6/YyvlnntTY9u3RS1LLdh88yhPHjnM84cljx9l98Gij2zfQS1LLtp69gfXr1rA24IR1a9h69oZGt2/qRpJadv6ZJ3HjFVvZffAoW8/e0GjaBgz0ktQJ5595UuMBfo6pG0kag3FW0YzKHr0kNWzcVTSjskcvSQ0bdxXNqAz0ktSwcVfRjMrUjSQ1bNxVNKMy0EvSGIyzimZUpm4kqXAGekkqnIFekgpnoJekwhnoJalwBnpJKpyBXpIKZ6CXpMIZ6CWpcAZ6SSqcgV6SCmegl6TCGeglqXAGekkqXK1AHxEXR8T+iDgQEdcOef3yiLizun0xIs6tu64kabyWDfQRsRa4DrgE2AK8ISK2LFjsm8DLM/Mc4F3AzhHWlSSNUZ0e/QXAgcw8mJlPADcB2+YvkJlfzMy5qc53A6fXXVeSNF51Av1pwP3zHh+unlvMbwL/Nuq6EbE9ImYjYvbIkSM1miVJzdp76DGu+8wB9h56bPmFe6TOVIIx5LkcumDEKxkE+peOum5m7qRK+czMzAxdRpLGZe+hx7j8+t08cew469et4cYrtnZmKsDVqtOjPwycMe/x6cCDCxeKiHOA64FtmXl0lHUlqW27Dx7liWPHOZ7w5LHj7D54dPmVeqJOoN8DbI6IsyJiPXApsGv+AhGxEfgE8MbMvHeUdSVpvrbSJ1vP3sD6dWtYG3DCujVsPXvDRN9/nJZN3WTmsYi4CrgVWAvckJn7IuLK6vUdwJ8AG4C/jwiAY5k5s9i6Y/osknquzfTJ+WeexI1XbGX3waNsPXtDMWkbqJejJzNvAW5Z8NyOefevAK6ou64kDTMsfTLJgHv+mSeN/f32Hnps4j8mtQK9JE3CXPrkyWPHi0ufQHtHLAZ6SZ1RcvoE2jtiMdBL6pRJpE/a0tYRi4FekiakrSMWA70kTVAbRywOUyxJhTPQS1LhDPSSVDgDvSQVzkAvaeqUOhzxYqy6kTRVmr46tY0hDUZloJc0VZq8OrUvY9ibupE0VZocjrgvY9jbo5c0VZq8OrUvg7BFZvdm7ZuZmcnZ2dm2myFJy+pKjj4i9mbmzLDX7NFL0ir0YRA2c/SSVDgDvSQVzkAvSYUz0EtqxbRdndomT8ZKmri+XGhUCnv0kiauLxcalcJAL2nimrw6VcszdSNNSFcurOmCtuZOnVYGemkCzEk/VR8uNCqFqRtpAsxJq00GemkC2s5JW8o43UzdSBPQZk66lLSR5zhWzkAvTUhbOemlJtroS/As5ceqLQZ6qXCLjZnep+DZ5KxQ08hALxVusbRRn4JnXyb46CoDvbQCfUl5zBmWNupT8LTufnWcYUoaUZ9SHstZ7Aerbz9kcoYpqVF9SnksZ1hPv6QfMg3UqqOPiIsjYn9EHIiIa4e8/uKIuD0iHo+Iaxa8dl9EfDUi7ogIu+nqvbZr4sfNi7vKs2yPPiLWAtcBrwEOA3siYldm3j1vsUeBtwGvW2Qzr8zMR1bbWKkLSs8X9yl3r3rqpG4uAA5k5kGAiLgJ2Ab8X6DPzIeBhyPiF8fSSqljSh6npfQfsmlUJ9CfBtw/7/Fh4MIR3iOB2yIigQ9l5s5hC0XEdmA7wMaNG0fYvKSmlfxDNo3q5OhjyHOjlOq8JDPPAy4B3hoRPz9soczcmZkzmTlzyimnjLB5SdJS6gT6w8AZ8x6fDjxY9w0y88Hq34eBmxmkgiRJE1In0O8BNkfEWRGxHrgU2FVn4xFxYkQ8a+4+cBFw10obK0ka3bI5+sw8FhFXAbcCa4EbMnNfRFxZvb4jIl4AzALPBo5HxNXAFuBk4OaImHuvj2bmJ8fzUSRJw9S6YCozbwFuWfDcjnn3/4tBSmeh7wDnrqaBkqTVceIRqWecRESjcggEqUccnkArYY9e6hGHJ9BKGOilHil9nB2Nh6kbqUe6OjyBwxp3m4Fe6pmuDU8wDecN+v5DZqCXtColjc8/TAk/ZOboJa1K6ecNSjgBbo9e0qp09bxBU0oYn985YyUVq6nceh9y9M4ZK2nqNJlb79oJ8FGZo5dUpBJy600x0EsqUukniUdh6kad0Yc8qPqj9JPEozDQqxNKqFVW9/Q9t94UUzfqhJXkU5scrtehf8fH77Z99ujVCaPWKjd5BODRxPj43XaDPXp1wlw+9e0X/UStYNBkRYXVGePjd9sN9ujVGaPkU5u8WrGEKx+7yu+2G7wyVmMxiQqaJt/Dip/x8budjKWujDXQq3HmZTVJ/pAMOASCJqr0YWvVHXYq6vFkrBrX9hWJlvNND0/21mOPXo1r84pEe3jTxZO99RjoNRZtXZFo2mi6OMxBPQZ6FcUe3vRxmIPlGehVFHt40lMZ6FWcvvXwLA/UuBnopRZ58liTYHml1CLLAzUJBnqpRW1fc6DpYOpGapEnjzUJBnqpZX07eaz+qZW6iYiLI2J/RByIiGuHvP7iiLg9Ih6PiGtGWVeSNF7LBvqIWAtcB1wCbAHeEBFbFiz2KPA24H0rWFeSNEZ1evQXAAcy82BmPgHcBGybv0BmPpyZe4AnR11Xkzetg35N6+eW6uToTwPun/f4MHBhze2vZl2NwbTWbU/r55agXo8+hjxXd7aS2utGxPaImI2I2SNHjtTcvEY1rXXb0/q5m+ZRUT/V6dEfBs6Y9/h04MGa26+9bmbuBHbCYIapmtvXiKZ10K9p/dxN8qiov+oE+j3A5og4C3gAuBS4rOb2V7OuxmBa67an9XM3ySGg+2vZQJ+ZxyLiKuBWYC1wQ2bui4grq9d3RMQLgFng2cDxiLga2JKZ3xm27rg+jOqZ1rrtaf3cTfGoqL+cHFxSbY602V1ODi6pER4V9ZODmklS4Qz0UkdZyqimmLqROshSRjXJHr3UQV7gpSYZ6KUGNZVuWcmEJKZ6tBhTNyMqvbys9M83Tk2mW0a9wMtUj5Yy1YF+1KBW+n+m0j/fuDV95egopYxetaqlTG3qZi6ovf+2/Vx+/e5ah7ul501L/3zj1ub8r849q6VMbY9+JT2g0i8BL/3zjVub4+k4lo+WMrVDIMz16OeCWt00Rek57NI/n1SqpYZAmNpADwY1Lc99RH3hWDeLcNyOyetT4PTktEox1YFek9W3wGkli0oxtVU306BrF9D0rarHShaVwh59obrYe+5bVY+VLCqFgb5QXUw79DFweh5HJTDQF6qrvWcDpzR5BvpC9bH3LGk8DPQFK6X33KeSTKmLDPTqtC6eVJb6xvJKdVrfSjKlLjLQq9OsZZdWz9SNOs2TytLqGejVeaWcVJbaYupGkgpnoG9J18ahkVQuUzct6GPJoLXsUn8Z6FvQxXFoltLHHyZJ/8/UTQv6VjJoLbvUb/boW9C3ksGuDpAmqZ6pnjNW9S2Wozd3L3WDc8Zq1YbVspu7l/rBHP0QTZY+llxGae5e6odaPfqIuBj4W2AtcH1m/uWC16N6/bXAD4Bfz8wvVa/dB3wX+CFwbLFDi65ospdaeo/X3L3UD8sG+ohYC1wHvAY4DOyJiF2Zefe8xS4BNle3C4EPVv/OeWVmPtJYq8eoydLHvpVRjqpvJ5WlaVWnR38BcCAzDwJExE3ANmB+oN8GfCQHZ3Z3R8RzIuLUzHyo8RaPWZO91Gno8ToOjdR9dQL9acD98x4f5kd764stcxrwEJDAbRGRwIcyc+ewN4mI7cB2gI0bN9Zq/Dg02Uu1xyupC+oE+hjy3MKazKWWeUlmPhgRzwM+FRFfy8zPPWXhwQ/AThiUV9Zo19g02Uu1xyupbXWqbg4DZ8x7fDrwYN1lMnPu34eBmxmkgrSIkqt0JLWjTqDfA2yOiLMiYj1wKbBrwTK7gDfFwFbg25n5UEScGBHPAoiIE4GLgLsabH9R5qp03n/bfi6/frfBXlIjlk3dZOaxiLgKuJVBeeUNmbkvIq6sXt8B3MKgtPIAg/LKt1SrPx+4eVB9yTrgo5n5ycY/RSFKr9KR1I5adfSZeQuDYD7/uR3z7ifw1iHrHQTOXWUbp8ZKq3QchkDSUhwCoUNWUqVT+kVZklbPQN8xo1bpmO6RtBzHuum5vo1tL2ny7NH3nBdlSVqOgb4AXpQlaSmmbiSpcAZ6SSqcgV6SCmegl6TCGeglqXAGekkqnIFekgpnoJekwhnoJalwBnpJKpyBXpIKZ6BviHO9SuoqBzVrgJN/SOoye/QNGDb5xzh41CBpJezRN2Clc72OwqMGSStloG/AJCb/cMpASStloG/IuCf/mMRRg6QyGeh7wikDJa2Ugb5HnDJQ0kpYdSNJhTPQS1LhpiLQW38uaZoVn6O3/lzStCu+Rz+pq1YlqauKD/Rz9edrA+vPJU2l4lM31p9LmnbFB3qw/lzSdCs+dSNJ085AL0mFqxXoI+LiiNgfEQci4tohr0dE/F31+p0RcV7ddSVJ47VsoI+ItcB1wCXAFuANEbFlwWKXAJur23bggyOs2xgvjJKkp6pzMvYC4EBmHgSIiJuAbcDd85bZBnwkMxPYHRHPiYhTgU011m2EF0ZJ0nB1UjenAffPe3y4eq7OMnXWBSAitkfEbETMHjlypEazfpQXRknScHUCfQx5LmsuU2fdwZOZOzNzJjNnTjnllBrN+lFeGCVJw9VJ3RwGzpj3+HTgwZrLrK+xbiO8MEqShqsT6PcAmyPiLOAB4FLgsgXL7AKuqnLwFwLfzsyHIuJIjXUb44VRkvRUywb6zDwWEVcBtwJrgRsyc19EXFm9vgO4BXgtcAD4AfCWpdYdyyeRJA0Vg0KZbpmZmcnZ2dm2myFJvRERezNzZthrXhkrSYUz0EtS4Qz0klQ4A70kFa6TJ2OrssxDbbcDOBl4pO1GDNHFdnWxTWC7RmW76utam87MzKFXm3Yy0HdFRMwudha7TV1sVxfbBLZrVLarvi62aTGmbiSpcAZ6SSqcgX5pO9tuwCK62K4utgls16hsV31dbNNQ5uglqXD26CWpcAZ6SSpcsYE+Is6IiM9ExD0RsS8ifrt6/rkR8amI+Hr170nV8xuq5b8XER9YZJu7IuKuJd7znIi4vXq/r0bE09tuV0ScEBEfrtpzT0S8Y9zfV0R8tpoQ/o7q9rxF3vMd1aTx+yPiF9puU0S8JiL2Vt/V3oh4VVe+q2rZjdU2rulKuya9z9f8O7axz6+PiJ0RcW9EfC0iXr/Iey65z49NZhZ5A04FzqvuPwu4l8EE5e8Frq2evxZ4T3X/ROClwJXAB4Zs71eBjwJ3LfJ+64A7gXOrxxuAtR1o12XATdX9ZwD3AZvG2S7gs8DMMn+fLcBXgKcBZwHfWPh9tdCmnwVeWN3/KeCBce9bddo1b9mPAx8DrulCu2hhn6/Zrjb2+XcC767urwFOXsk+P65bsT36zHwoM79U3f8ucA+D+Wq3AR+uFvsw8Lpqme9n5ueB/1m4rYh4JvB24N1LvOVFwJ2Z+ZVqe0cz84cdaFcCJ0bEOuDHgCeA74yzXTVtY/Cf8fHM/CaDuQwuaLNNmfnlzJybAW0f8PSIeNqQ5Sb9XRERrwMOVu1arP2TbtfE9/ma2tjnfwP4i2q545k57IrZZff5cSk20M8XEZsY9Nb+A3h+Zj4Egz80sOih8jzvAt7PYFKVxbwIyIi4NSK+FBG/15F2/RPwfeAh4FvA+zLz0TG3C+Afq0PrP46IYXMH1544foJtmu/1wJcz8/GlFppEuyLiROD3GfQaa5nQ99XGPl+nXRPd5yPiOdXdd1Xfw8ci4vlDFh1pn29S8YG+6vV+HLg6M5/yq15j/Z8Bfjwzb15m0XUMDusur/79lYh4dQfadQHwQ+CFDA4Xfycizh5XuyqXZ+ZPAy+rbm8c9lZDnhta6zvBNs29308C7wF+a6kNTrBd7wT+OjO/V2eDE2zXRPf5Edo16X1+HYP5sL+QmecBtwPvG/ZWQ56bSH170YE+Ik5g8Ae8MTM/UT393xFxavX6qcDDy2zm54DzI+I+4PPAiyLis0OWOwz8e2Y+kpk/YDC94nkdaNdlwCcz88nMfBj4AjB0fI6G2kVmPlD9+10G5w+GHZ7WmXR+0m0iIk4HbgbelJnfWGx7E27XhcB7q7/11cAfxGCKzrbbNel9vm67Jr3PH2VwVD3X6foYw7+HWvv8OBQb6KtDun8A7snMv5r30i7gzdX9NwP/vNR2MvODmfnCzNzEoNdyb2a+YsiitwLnRMQzqtzgy4G7O9CubwGvioETga3A18bVrohYFxEnV/dPAH4JGFYRtAu4NCKeFoPJ4zcD/9lmm6pD8H8F3pGZX1hiexNtV2a+LDM3VX/rvwH+PDOfUoHVwt9wovv8CO2a6D6fmQn8C/CK6qlXM+R7oMY+PzY5gTO+bdwYBL9kUBVwR3V7LYPKgE8DX6/+fe68de4DHgW+x+DXd8uCbW5iXnUL8MvAn817/GsMTpbdBby3C+0Cnsmgh7GPwc73u+NsF4PKhL3VdvYBf0tVWTDk+/pDBpUH+4FL2m4T8EcMcrt3zLs9r+12LXjvP2Xxqps2/oYT2+dH+DtOdJ+vnj8T+Fy1rU8DG1eyz4/r5hAIklS4YlM3kqQBA70kFc5AL0mFM9BLUuEM9JJUOAO9JBXOQC9Jhftf0tLRSoQvMakAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "import numpy as np\n",
    "from scipy import linalg, stats\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "## simulate displacement time-series and dates\n",
    "vel_sim = 0.1 #m/year\n",
    "num_date = 50\n",
    "t_year = np.arange(num_date) * 12 / 365.25 + 2014.5\n",
    "\n",
    "np.random.seed(65535)\n",
    "dis_sim = (t_year - t_year[0]) * vel_sim + np.random.rand(num_date) * 0.1\n",
    "dis_sim -= dis_sim[0]\n",
    "\n",
    "print('Displacement [m]:', dis_sim)\n",
    "print('Dates [year]:', t_year)\n",
    "print('True velocity [m/year]:', vel_sim)\n",
    "\n",
    "# plot\n",
    "plt.figure()\n",
    "plt.plot(t_year, dis_sim, '.')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Estimation using scipy.linalg with velocity [m/year]: 0.10, velocity std [m/year]: 0.01\n",
      "Estimation using scipy.stats  with velocity [m/year]: 0.10, velocity std [m/year]: 0.01\n",
      "LinregressResult(slope=0.09950559842405902, intercept=-200.42318059638643, rvalue=0.8588459418905924, pvalue=1.5017092073839943e-15, stderr=0.008566038308837858)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de3xU1bn/8c+TEAgkXBQCgoAgoMgtAcLNSim1CCIFKVqsqKhYRGtbe079Fdvf8bT2nF89tcfWO6LoQYvFSwviUeultmrlGjQiN5VQwABKuCfck6zfH3ugAWYnk2SuO9/36zUvMjPP3vuZzPBkzdprr2XOOUREJLjSEp2AiIjElgq9iEjAqdCLiAScCr2ISMCp0IuIBFyjRCcQTps2bVyXLl0SnYaISMpYuXLlTudcTrjnkrLQd+nShYKCgkSnISKSMsxss99z6roREQk4FXoRkYBToRcRCbik7KMP59ixYxQXF3P48OFEpyJxlpmZSceOHcnIyEh0KiIpKWUKfXFxMc2bN6dLly6YWaLTkThxzrFr1y6Ki4vp2rVrotMRSUkp03Vz+PBhWrdurSLfwJgZrVu31jc5kXpImRY9oCLfQOl9l6Dbsusg0+auYGPJAc7NyWLO1EF0bt0savtPmRa9iEhQTZu7gqKSMiqco6ikjGlzV0R1/yr0tZCenk5eXt6J2z333OMbu3DhQtauXXvi/l133cVbb71V7xz27t3LI488Uuvtfv7zn/Ob3/ym3sev6362bdvGFVdcAUBhYSGvvvpqvXMRCYqNJQeoDC0NUum8+9GUUl03ida0aVMKCwsjil24cCHjxo2jV69eANx9991RyeF4ob/11lujsr946dChAy+++CLgFfqCggLGjh2b4KxEksO5OVkUlZRR6SDNvPvRlHot+p+3jO2tDmbOnEmvXr3o168fP/7xj1m8eDGLFi3ijjvuIC8vj6KiIq6//voTha5Lly789Kc/ZdiwYeTn5/PBBx8wevRounXrxqxZswAoKyvj4osvZsCAAfTt25eXXnrpxLGKiorIy8vjjjvuAODee+9l0KBB9OvXj3//938/kdd//ud/cv755/ONb3yDTz755LS89+3bR5cuXaisrATg4MGDdOrUiWPHjlFUVMSYMWMYOHAgw4cPZ/369adtX1hYyNChQ+nXrx8TJ05kz549AGzYsIFvfOMb5ObmMmDAAIqKiti0aRN9+vTh6NGj3HXXXTz33HPk5eXx3HPP0aNHD0pKSgCorKyke/fu7Ny5s07vhUgqmjN1EN1yskk3o1tONnOmDorq/tWir4VDhw6Rl5d34v6dd97JqFGjWLBgAevXr8fM2Lt3L61atWL8+PGMGzfuRHfFqTp16sSSJUv40Y9+xPXXX8/777/P4cOH6d27NzNmzCAzM5MFCxbQokULdu7cydChQxk/fjz33HMPq1evPvHN4o033uCzzz5j+fLlOOcYP3487777LllZWcyfP58PP/yQ8vJyBgwYwMCBA0/KoWXLluTm5vLOO+8wcuRIXn75ZUaPHk1GRgbTp09n1qxZ9OjRg2XLlnHrrbfy9ttvn7T9ddddx4MPPsiIESO46667+MUvfsHvfvc7pkyZwsyZM5k4cSKHDx+msrKSHTt2ANC4cWPuvvtuCgoKeOihhwBYv3498+bN4/bbb+ett94iNzeXNm3aRO19E0l2nVs3481/GRGz/avQ10K4rpvy8nIyMzO56aabuOyyyxg3blxE+xo/fjwAffv2paysjObNm9O8eXMyMzPZu3cvWVlZ/PSnP+Xdd98lLS2NrVu38uWXX562nzfeeIM33niD/v37A943gc8++4zS0lImTpxIs2bNTjreqSZPnsxzzz3HyJEjmT9/PrfeeitlZWUsXryYK6+88kTckSNHTtpu37597N27lxEjvA/n1KlTufLKKyktLWXr1q1MnDgR8C52qsmNN97IhAkTuP3223nyySe54YYbatxGJJnFehRNbaVe102SadSoEcuXL2fSpEksXLiQMWPGRLRdkyZNAEhLSzvx8/H75eXlzJs3j5KSElauXElhYSHt2rULO5bcOcedd95JYWEhhYWFbNiwgWnTpgGRDUscP348r732Grt372blypV8/etfp7KyklatWp3YZ2FhIevWrYvoddVlsflOnTrRrl073n77bZYtW8all15a632IJJNYj6KprdRr0f98X6IzOElZWRkHDx5k7NixDB06lO7duwPQvHlzSktL67zfffv20bZtWzIyMvjrX//K5s2bw+539OjR/Nu//RtTpkwhOzubrVu3kpGRwVe/+lWuv/56Zs6cSXl5OS+//DI333zzacfJzs5m8ODB/PCHP2TcuHGkp6fTokULunbtygsvvMCVV16Jc45Vq1aRm5t7YruWLVtyxhln8N577zF8+HCeeeYZRowYQYsWLejYsSMLFy7k8ssv58iRI1RUVJx0zHC/m5tuuolrrrmGa6+9lvT09Dr/3kSSQaxH0dSWWvS1cLyP/vht5syZlJaWMm7cOPr168eIESP47W9/C8BVV13FvffeS//+/SkqKqr1saZMmUJBQQH5+fnMmzePnj17AtC6dWu+8pWv0KdPH+644w4uueQSrr76aoYNG0bfvn254oorKC0tZcCAAUyePJm8vDwmTZrE8OHDfY81efJkfv/73zN58uQTj82bN485c+aQm5tL7969T5wMrmru3Lnccccd9OvXj8LCQu666y4AnnnmGR544AH69evHhRdeyBdffHHSdiNHjmTt2rUnTsaC982irKxM3TYSCOfmZJEW+kIdi1E0tWV1+aoda/n5+e7UhUfWrVvHBRdckKCMJNYKCgr40Y9+xHvvvRf2eb3/kkoS0UdvZiudc/nhnku9rhsJnHvuuYdHH32UefPmJToVkaiI9Sia2lLXjSTczJkz2bx5MxdddFGiUxEJpJQq9MnYzSSxp/ddpH5SptBnZmaya9cu/advYI7PRx/JeHwRCS9l+ug7duxIcXHxiUvlpeE4vsKUiNRNyhT6jIwMrTAkIlIHKdN1IyIidaNCLyIScCr0IiIBp0IvIhJwKvQiIgGnQi8iEnAq9CIiAadCLyIScCr0IiIBF1GhN7MxZvaJmW0ws5lhnp9iZqtCt8VmlhvptiIiEls1FnozSwceBi4FegHfMbNep4T9AxjhnOsH/BKYXYttRUQkhiJp0Q8GNjjnNjrnjgLzgQlVA5xzi51ze0J3lwIdI91WRERiK5JJzc4GPq9yvxgYUk38NOC12m5rZtOB6QCdO3eOIC0RkehJxPJ/8RJJi97CPBZ2UngzG4lX6H9S222dc7Odc/nOufycnJwI0hIRiZ5pc1dQVFJGhXMUlZQxbe6KRKcUNZG06IuBTlXudwS2nRpkZv2AJ4BLnXO7arOtiEiibSw5QGWoGVrpvPtBEUmLfgXQw8y6mllj4CpgUdUAM+sM/Am41jn3aW22FRE5bsuug4y67x263fkqo+57hy27Dsbt2OfmZJEW6oNIM+9+UNRY6J1z5cBtwOvAOuB559waM5thZjNCYXcBrYFHzKzQzAqq2zYGr0NEAiCR3Sdzpg6iW0426WZ0y8lmztRBcTt2rFkyrsGan5/vCgoKEp2GiMRZtztfpaJKTUo3o+hXYxOYUXRVe8K34hisfQna9ID2udXvKAwzW+mcyw/3nK6MFZGkEeTuE/D5xnJ4Hyx+EO7Pgz9Og3fvjfpxVehFJGkEufsETj7h296VcNXuWXBfb3jj/8L+Yu+Jdf8LuzdG9bgpszi4iARf59bNePNfRiQ6jZg5NyeLrJ0fMS39FS5NW04jq4Sjp0Y5WPoojI1ey16FXkQk1ior4JPXeDnrfjL3L/ePa5QJeVfD4OlRPbwKvYhIrBw9CIXzYOkjsHsjmX5xWTlecc+fBlmto56GCr2ISLSVfgkrHocVT8ChPf5xOT1h2Peg77chw/fPQL2p0IuIRMuXa2HJw/Dx81BxWuf7P3UdARf+ALpfDBZuppjoUqEXEakP52DjX2HxQ1D0F/+4tAzoe4XXgj+rb/zyQ4VeRKRuyo/C6he9FvyXq/3jMltC/o1eH3yLDvHLrwoVehFpUOo9HfHB3bDyKVg2G8q+8I9rdY7Xes+bAk2y6594PajQi0iDcvzq1ErHiatTIxq7v3ujN779w9/DsWomW+s4GC68DXqOg7T06CVeDyr0ItKg1Go6Yufg82XeFAXrX8FnOQ0qnPFG5SBea34FD9x0c/STricVehFpUM7NyTrRovedT6eiHNa/7J1g3eo/weIB14TnK77GkxVj+Ny1I32P8UAMc68rFXoRaVDmTB10Wh/9CUdKva6ZpY/A3i3+O2neHobczDXLzuOjnVT/RyMJqNCLSIMSdj6dfVth+WNQ8D9wZJ//xmf1hWHfh94ToVFj7r/g9BO7yUiFXkQaru0fed0za/4EleX+cT0ugWG3QdevnnSBU6pMwqZCLyINS2UlbHjTO8G66T3/uPQmkDsZhn4P2vaMX34xoEIvIg3DsUOw6jnvAqedn/rHNWsNg27ybtlt45dfDKnQi0iwlZV4k4uteAIO7vSPa93du8Ap9zuQ0TR++cWBCr2IxF29r06NRMmnsOQh+Gg+VBzxjzvnIu8Cpx6jIS2Yi+6p0ItI3NX56tSaOOf1uy9+CD573T/O0r2RMxfeBh361/+4SU6FXkTirlZXp0ai4hisWeCdYP1ilX9ckxYw4DoYMgNadarfMVOICr2IxF1EV6dG4tBe+GAuLHsM9m/1j2vZCYbeAv2vhcwWdTtWClOhF4mDuPRJp5Bqr06NxJ7NoQnGnoGjZf5xHQZ4J1h7XQ7pDbfcmXPhJ+lJpPz8fFdQ4D+/hEiqGXXfOye1YLvlZKfEhTZJp7jA655ZtwhcpU+Qwfljvf73zsPisoJTMjCzlc65/HDPNdw/cSJxFPU+6YaksgI+edU7wfr5Uv+4Rk2h/xQYeiu07ha//FKACr1IHEStT7oOUrbb6OgBKHzWm2Bs90b/uKy2MGQ65E+DZmfGL78UokIvEgf17pOuh5gNZYyV0i9g+WxYMQcO7/WPa9vL63/veyU0ahK//FKQCr1IHCRy8iu/bqOka+l/ucabnmDV81B5zDfsvcq+vJo9iV/dcnuD6X+vLxV6kYDz6zZKipa+c1D0F6//feNffcOOunReqvgKT1SM5RPXmfQ9xq9U5COmQi8ScH7dRgk9QVx+BD5+wWvB71jrH5fZCgZN49oP+7JiV+OkX+AjWanQi9RS0nV51MCv2yghJ4gP7oaCObD8cSj70j/ujK5e/3ve1dA4i3vzUmOBj2SlcfQitRSUMfF+f7Bi8odsV5E3eubDeVB+yD+u01Bv/Pv5YyEtvX7HbGA0jl4kioIyJt6vpR+1vnvnYMtSbwbJ9a8APo1KS4NeE7wVnDqGrVNSTxEVejMbA9wPpANPOOfuOeX5nsBTwADgZ86531R5bhNQClQA5X5/cURSRSLHxMdDvf+QVZTDupe8E6zbPvCPa5z9zwnGzjin7glLjWos9GaWDjwMjAKKgRVmtsg5V/UMym7gB8DlPrsZ6ZyrZsZ/kdSRyDHx8VDnP2SH93tzzyydBfu2+Me1OBuG3AwDpkLTVtFJWqoVSYt+MLDBObcRwMzmAxOAE4XeObcD2GFml8UkS5EkkioLQtdVrf+Q7SuGZbNg5Vw4st8/7qx+cOH3vXng0zOim7RUK5JCfzbweZX7xcCQWhzDAW+YmQMec87NDhdkZtOB6QCdO3euxe5FJJoi/kO2rdDrf1+zACrL/eN6jPZOsHYZrgucEiSSQh/unanNUJ2vOOe2mVlb4E0zW++ce/e0HXp/AGaDN+qmFvsXkXiprITP3vAK/Kb3/OPSm0DuVd4QyZzz45efhBVJoS8Gqi7F0hHYFukBnHPbQv/uMLMFeF1BpxV6EUlixw7BR3+AJY/Ars/845q1hkHfhUE3QXZO/PKTakVS6FcAPcysK7AVuAq4OpKdm1kWkOacKw39fAlwd12TFZE4KyuBFY/Diifg4C7/uDbnea33fpMho2n88pOI1FjonXPlZnYb8Dre8MonnXNrzGxG6PlZZnYWUAC0ACrN7HagF9AGWGBev1wj4Fnn3J9j81JEJGp2rPe6Z1Y9DxVH/OO6DPdOsHYfBWlp8ctPaiWicfTOuVeBV095bFaVn7/A69I51X4gtz4JikicOAf/eMcb/77hTf+4tEbQZ5LXgm+v/96pQFfGiqSQmExPUH4U1vzJa8F/8bF/XJOWkH89DL4ZWp5dv2NKXKnQi6SQqE4tfGgvrHwKlj0Gpdv941p19pbn638NNGlet2NJQqnQi6SQqMyzs2cTLH0UPngGjlWz/dn53vj3nt+EdJWKVKZ3TySF1Guenc9XwJIHYd3L4Cp9ggx6XuadYO00RBc4BYQKvUgKqfX0BJUV3syRSx6Cz5f5x2U0g7wpMPQWaN2tVjml2vz8DZHmoxcJoiNlUDjPmwN+zyb/uOx2MHg65N8Izc6s06GCMj+/n1T5Q6b56EUaiv3bYfljUPAUHN7rH9e2t9f/3mcSNGpSr0MGZX5+P0mxtm49qdCLBMEXH3vrr378IlQe84/rdrFX4M8dGbX+d83Pn/xU6EVSlXOw4S1Y/KB3oZOf9MbQ99veBU7tekU9Dc3Pn/zURy+Sao4dho+f91rwJev945qe4U0uNui70Lxd/PJLEtHqWw9CH70KvUiqOLALCubA8sfhwA7/uDO7wbBbIfdqaJx8BSlegn6S+FQ6GSuSynZugKUPQ+EfoPyQf1znC73+9/Mu1QRjBKNvPVpU6EWSkXOwebE3/v2T1/Bd68fSodcEr8CfPTCuKSa7IPStR4sKvSSFVOkHjbmKY7D2Ja/Ab/vQP65xcxhwHQyd4c1FI6cJ+kni2lAfvSSFhtafeprD++CDp70JxvZ97h/XoqNX3AdcB5kt45efJD310UvSq21/ajS/AST028TeLbB0llfkj5b6x7XP8+af6TUB0jPik1sU6JtactAZG0kK5+ZkkRa6fieS/tTjVytWOHfiasW6iua+IrZ1JbxwA9yf551o9Svy510K178C0/8Gfa9IqSIPCfrdymnUopekUNv+1GiOqIjb6IzKSvj0NW8Fpy2L/eMaZULud7wLnNr0iE0ucaKRL8lBhV6SQufWzWrVJx/NERUxH51x9CB89CwseQR2F/nHZeV4FzcNmgZZbaKbQ4Jo5EtyUNeNRN2WXQcZdd87dLvzVUbd9w5bdh2M+jHmTB1Et5xs0s3olpNdrxEV0dzXSUq/hLf/A37bG175V/8i3+Z8+OYDcPtq+NpPAlPkIYa/W6kVjbqRqGvwI2h2rPOGR656HiqO+sd1HeGdYO12sS5wqiOd7P0njbqRuGqQ/bLOwca/eQV+w1v+cWmNoM8VXv97+35xSy+ogjCFcDyo0EvUJbJfNu4tvPKjsPpFb4KxL1f7x2W2hIE3wJCboUWH2OXTwDTIRkUd6PuiRF0i+2XjNpzv4G5477/hd31h4S3+Rb7VOTDmv+BHa2HUL1Tko6y2w3IbKrXoJepqO4ImmmLewtu9EZY+Ch/+Ho5Vc5K54yAYdhtc8E1IS49uDnKCpjmIjAq9BErMuo22LIMlD8K6/8V/grE06DnOO8HaaXB0jivVSmSjIpWo0EugRLWFV1EO61/2+t+Lq+kCysiC/tfA0FvgzK51P55IjKjQS6BEpYV3pNTrmln6iDcXjZ/m7b2TqwOv91ZzqgMND5R4UKEXOW7fVlj+GBT8DxzZ5x/Xrq83/3vvb0GjxvU6pIYHSjyo0ItsX+WNf1/9R6gs94/rcYl3grXrV8EsKofW8ECJBxV6aZgqK70Lm5Y8CP941z8uvQn0+7ZX4Nv2jHoamgtG4kGFXhqWY4dh1XxvgrGdn/jHNWsNg27ybtltY5aOhgdKPKjQS8NwYCeseAKWPw4Hd/rHte7uTU+Q+x3IaBrztDQ8UOIhokJvZmOA+4F04Ann3D2nPN8TeAoYAPzMOfebSLcViamST72FPT6aD+WH/ePOucg7wdpjtCYYk8CpsdCbWTrwMDAKKAZWmNki59zaKmG7gR8Al9dhW5Hocg42/d07wfrpn/3jLB16T/QKfIf+8ctPJM4iadEPBjY45zYCmNl8YAJwolg753YAO8zsstpuK/EV6HHbFcdgzULvBOv2j3zDDtCMiv5TaTHie9CqUxwTFEmMSL6jng1UXZa+OPRYJOqzrcRAINfwPLwP3n8A7s+FP93kW+S3ujb88tg1DDvyAJOKLlWRlwYjkhZ9uAHDka5WEvG2ZjYdmA7QuXPnCHcvtRWocdt7NsOyWfDB03C0zD+uwwB+sPkiXqkYRAXeBGMHUvl1J0igvw0GXCQt+mKgatOnI7Atwv1HvK1zbrZzLt85l5+TkxPh7qW2AjGta/FKeOF6eKC/N01B2CJvcP5lcMNr8N23Wdf6GzjzinzKvu4EC+S3wQYikkK/AuhhZl3NrDFwFbAowv3XZ1uJgZRdw7Oywps58skx8MTXYc0CcBWnxzVqCvnT4LYC+M6zcM6FYJa6rzuJBOrbYANTY9eNc67czG4DXscbIvmkc26Nmc0IPT/LzM4CCoAWQKWZ3Q70cs7tD7dtrF6M1Czlxm0fPQCFz3ot990b/eOy2sLg6ZB/I2S1Pu3plHvdSUhX8aYuLQ4uyan0C1g+GwqehEN7/ONyLvAucOp7JWRkxi+/Bkh99MlNi4NL6vhyjTf/+8cvQMVR/7hzR3rj37tdHLUJxqR6+laUulToJfGcg6K3vQucit72j0vL8Fruw74HZ/WJX34iKU6FXhKn/IjXcl/yMOyo5hq6zFZe3/vg6dCiffzySyB1k0g0qdBL/B3cDQVzvAnGyr70jzujq9d6z7saGjesE39akESiSYVe4mdXkTd6pvBZOHbQP67TUK///fyxkJYev/ySiIYySjSp0EtsOQdblnr97+tfwfeiakuDC8bDhd+HjmEHDiS9aHa31HYoo7p6pDoaXlkLQf/PFNXXV1EO6xZ5BX7rSv+4xtnQ/1oYOgPO6FK3YyWJUfe9c1Jx7paTXefultq+F9E8tqQmDa8Moy5FLej9plF5fUdK4YNnYNmjsHeLf1zzDjDkZhh4PTRtVa+8k0U0u1tqO5RRXT1SnQZb6OtS1IL+n6ler29fsTfB2Mq5cGS/f9xZfWHY97154Bs1rl/CSSaRV47qqlWpToNdSqcuRS0QE4JVo06vb1sh/PEmb4rgxQ/6F/keo2Hqy3Dze5A7OXBFHhI7j5Dm8pHqNNg++rr0aaqPPqSyEj57HRY/BJv/7r/D9CaQe5U3RDLn/NglLiLV9tE32EIf9KIdE8cOwUd/gCWPwK7P/OOatYHB3/VmkcxO3Smn9RmRVKJCL/VTVgIrHocVT8DBXf5xbc7zWu/9JkNG09OeTrXCqZEskko06kbqZsd6WPowfPQcVBzxj+sy3Bv/3n0UpPmf9km1UUtBP/kuDYcKfUDVufXsHPzjHa//fcOb/nFpjaD3t7wWfIe8iHJKtcKpkSwSFA121E3Q1XrZt/KjXsv9seHw9AT/It+kJVz4A/jhKpj0eMRFHlJv1JJGskhQqEUfUBG3ng/thZVPwbLZUFrNUsCtOsOQW2DAtdCkeZ1ymjN10GnfMpKZ5l+XoFChD6gaux32bIKlj3pXsR6rpgvl7HxvgrGe34T0+n1cVDhFEkOFPqB8W8+fr4AlD8K6l8FV+mxt0PMy7wRrpyFawUkkxanQB9RJrefKCm/myIUPwefL/DfKaAZ5U2DoLdC6W3wSrUGqDckUSUYq9EF29AB8OM+bA37PP/zjss+CIdNh4A3Q7Mz45ReBVBuSKZKMVOiDaP92WD4bCp6Ew3v949r29vrf+0yCRk3il18tpNqQTJFkpEIfJF987K2/+vGLUHnMP67bxV6BP3dk0ve/ayy7SP2p0Kc652DDW94CHxv/5h+X3hj6ftu7wKldr7ilV1+pNiRTJBmp0KeqY4fh4+e9FnzJev+4pmd4k4sNng7N28UvvyjRkEyR+lOhTzUHdkHBHK8P/kCJf9yZ58LQWyHvamis7g6RhkyFPgHqNGRw5wZvgrHCP0D5If+4zhd6/e/njYG09OgmLiIpSYU+ASIeMugcbF7s9b9/8hrgM6W0pUOv8d4SfR0HRj1fjWUXSW0q9AlQ45DBinJYu9Ar8Ns+9N9R42wYMNVbZPuMc2KWr8ayi6Q2FfoE8B0yeHg/fPC0t8j2vs/9d9DibBgyAwZOhcyWMc9XY9lFUpsKfQKcOmTwf77VHl7/GaycC0dL/Tdsn+t1z/S+HNIz4pavxrKLpDYtJZhIWz/wumfWLARX4R933hgYdht0uSghFzj59dGr714keWjN2GRSWQmf/tkr8Jvf949rlAm53/GGSOacF7/8akFrqookD60ZWwvRbKVW3dcFbRrx9MAizlz1OOwu8t8oKwcGfRcGTYOsNnV8FfGhvnuR1BDRUoJmNsbMPjGzDWY2M8zzZmYPhJ5fZWYDqjy3ycw+NrNCM0v6Znqtl+CrYV/7Srbyw/TneHr/DZz5t5n+Rb7N+fDNB+D21fC1nyR9kYfUWxpQpKGqsUVvZunAw8AooBhYYWaLnHNrq4RdCvQI3YYAj4b+PW6kc25n1LKOoai1Unes47t77mNC47/TxMr947qO8Bb46HYxpKXWEr6ah0YkNUTSdTMY2OCc2whgZvOBCUDVQj8BeNp5Hf5LzayVmbV3zm2PesYxVq8RJs55E4steQg2vMW3/S5MTWsEfa7wJhhr3y8aaSeE5qERSQ2RFPqzgaqDuos5ubXuF3M2sB3vcs43zMwBjznnZoc7iJlNB6YDdO7cOaLkY6FOrdTyo7D6j94EY19+7BtW2aQlafk3eBc4tegQxaxFRPxFUujDjec7dahOdTFfcc5tM7O2wJtmtt459+5pwd4fgNngjbqJIK+YqFUr9dAeKHjKm2CstJovL63OgaG3ktb/GmiSHZ1ERUQiFEmhLwY6VbnfEdgWaYxz7vi/O8xsAV5X0GmFPqXs3ghLH4UPfw/HDvrHdRzkjX+/4JsRTTCmcekiEguRnJvp0UkAAAnbSURBVP1bAfQws65m1hi4Clh0Sswi4LrQ6JuhwD7n3HYzyzKz5gBmlgVcAqyOYv7x9flyeO4aeHCg14oPV+QtDS4YD9PehJve8q5ijXAWyWiO+BEROa7GFr1zrtzMbgNeB9KBJ51za8xsRuj5WcCrwFhgA3AQuCG0eTtggXlXczYCnnXO/TnqryKWKitg3cveCdbiagpvRjPofw0MvcWbC74ONC5dRGIhogumnHOv4hXzqo/NqvKzA74XZruNQG49c0yMI6Xw4TxY+gjs3ewf17y9t3pT/g3eak71UJcRP+ruEZGa6MrYU+3fBsseg5VPweF9/nHt+nj9730mQaPGUTl0XUb8aAphEamJCv1x21d53TOr/wiV1Vzg1H2UN/793K9FfYKxuoxLV3ePiNSkYRf6ykrY8BYseRD+Uc1AoPTG0G+y14Jv2zN++UVAUwiLSE0aZqE/dhhWPef1v5es949reiYMugkGfxey28Yvv1rQNAQiUpOGVegP7IQVc2DF43CgxD+udXdveuDc70Dj5D6xqWkIRKQmDaPQ7/zM63//aD6UH/aPO+ciuPA26DE65SYYExHxE9xC7xxs+rtX4D+tZui+pXsXNQ27Dc4e4B8nIpKiglfoK455S/MteRC2f+Qf17i5t7j2kBnQqpN/nIhIigtOoT+8z1tce9ks2L/VP65lJ6+4D7gOMlvELz8RkQQJRqH/+2/h3f+Go6X+MR36e90zvS6H9GC8bBGRSASj4lm6T5E3OH+sd4HTORdG/QInEZFUEIxCP3AqvPPrfxb7Rk0h72pviGSb7onNTUQkwYJR6DNben3uH78QmmDsRshqneisRESSQjAKPcDXfgIX3wUZmYnOREQkqQSn0Ge2THQGIiJJSZd/iogEXHBa9AmkxT9EJJmpRR8FWutVRJKZWvRREI/FP/StQUTqSi36KDg3J4u00LVYsVr8Q98aRKSuVOijYM7UQXTLySbdjG452TFZ/ENLBopIXanrJgrisfiHlgwUkbpSiz5FxONbg4gEk1r0KUJLBopIXalFLyIScCr0IiIBF/iuG40/F5GGLvAteo0/F5GGLvCFXuPPRaShC3yhj8dVqyIiySzwhV7jz0WkoQv8yViNPxeRhi7wLXoRkYZOhV5EJOAiKvRmNsbMPjGzDWY2M8zzZmYPhJ5fZWYDIt1WRERiq8Y+ejNLBx4GRgHFwAozW+ScW1sl7FKgR+g2BHgUGBLhtlGhC6NERMKLpEU/GNjgnNvonDsKzAcmnBIzAXjaeZYCrcysfYTbRoUujBIRCS+SQn828HmV+8WhxyKJiWRbAMxsupkVmFlBSUlJBGmdTBdGiYiEF0mhtzCPuQhjItnWe9C52c65fOdcfk5OTgRpnUwXRomIhBdJoS8GOlW53xHYFmFMJNtGhS6MEhEJL5ILplYAPcysK7AVuAq4+pSYRcBtZjYf72TsPufcdjMriWDbqNCFUSIi4dVY6J1z5WZ2G/A6kA486ZxbY2YzQs/PAl4FxgIbgIPADdVtG5NXIiIiYZlzYbvMEyo/P98VFBQkOg0RkZRhZiudc/nhntOVsSIiAadCLyIScCr0IiIBp0IvIhJwSXkyNjQsc3Oi8wDaADsTnUQYyZhXMuYEyqu2lFfkki2nc5xzYa82TcpCnyzMrMDvLHYiJWNeyZgTKK/aUl6RS8ac/KjrRkQk4FToRUQCToW+erMTnYCPZMwrGXMC5VVbyityyZhTWOqjFxEJOLXoRUQCToVeRCTgAlvozayTmf3VzNaZ2Roz+2Ho8TPN7E0z+yz07xmhx1uH4svM7CGffS4ys9XVHLOfmS0JHe9jM8tMdF5mlmFmc0P5rDOzO2P9+zKzv4UWhC8M3dr6HPPO0KLxn5jZ6ETnZGajzGxl6He10sy+niy/q1Bs59A+fpwsecX7Mx/h+5iIz3xjM5ttZp+a2Xozm+RzzGo/8zHjnAvkDWgPDAj93Bz4FOgF/BqYGXp8JvBfoZ+zgIuAGcBDYfb3LeBZYLXP8RoBq4Dc0P3WQHoS5HU1MD/0czNgE9AllnkBfwPya3h/egEfAU2ArkDRqb+vBOTUH+gQ+rkPsDXWn61I8qoS+0fgBeDHyZAXCfjMR5hXIj7zvwD+I/RzGtCmLp/5WN0C26J3zm13zn0Q+rkUWIe3Xu0EYG4obC5weSjmgHPu78DhU/dlZtnAvwD/Uc0hLwFWOec+Cu1vl3OuIgnyckCWmTUCmgJHgf2xzCtCE/D+Mx5xzv0Dby2DwYnMyTn3oXPu+Apoa4BMM2sSJi7evyvM7HJgYygvv/zjnVfcP/MRSsRn/kbgV6G4SudcuCtma/zMx0pgC31VZtYFr7W2DGjnnNsO3hsN+H5VruKXwH/jLari5zzAmdnrZvaBmf2fJMnrReAAsB3YAvzGObc7xnkBPBX6av1vZhZu7eCIF46PY05VTQI+dM4dqS4oHnmZWRbwE7xWY0Ti9PtKxGc+krzi+pk3s1ahH38Z+j28YGbtwoTW6jMfTYEv9KFW7x+B251zp/1Vj2D7PKC7c25BDaGN8L7WTQn9O9HMLk6CvAYDFUAHvK+L/2pm58Yqr5Apzrm+wPDQ7dpwhwrzWNixvnHM6fjxegP/Bdxc3Q7jmNcvgN8658oi2WEc84rrZ74WecX7M98Ibz3s951zA4AlwG/CHSrMY3EZ3x7oQm9mGXhv4Dzn3J9CD39pZu1Dz7cHdtSwm2HAQDPbBPwdOM/M/hYmrhh4xzm30zl3EG95xQFJkNfVwJ+dc8ecczuA94Gw83NEKS+cc1tD/5binT8I9/U0ooXj45wTZtYRWABc55wr8ttfnPMaAvw69F7fDvzUvCU6E51XvD/zkeYV78/8Lrxv1ccbXS8Q/vcQ0Wc+FgJb6ENf6eYA65xz91V5ahEwNfTzVOCl6vbjnHvUOdfBOdcFr9XyqXPua2FCXwf6mVmzUN/gCGBtEuS1Bfi6ebKAocD6WOVlZo3MrE3o5wxgHBBuRNAi4Coza2Le4vE9gOWJzCn0FfwV4E7n3PvV7C+ueTnnhjvnuoTe698B/885d9oIrAS8h3H9zNcir7h+5p1zDngZ+FrooYsJ83sggs98zLg4nPFNxA2v+Dm8UQGFodtYvJEBfwE+C/17ZpVtNgG7gTK8v769TtlnF6qMbgHGA3dXuX8N3smy1cCvkyEvIBuvhbEG78N3RyzzwhuZsDK0nzXA/YRGFoT5ff0Mb+TBJ8Clic4J+L94fbuFVW5tE53XKcf+Of6jbhLxHsbtM1+L9zGun/nQ4+cA74b29Regc10+87G6aQoEEZGAC2zXjYiIeFToRUQCToVeRCTgVOhFRAJOhV5EJOBU6EVEAk6FXkQk4P4/Hp1wmQYvBfAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## least squares solution\n",
    "# option 1 - scipy.linalg module\n",
    "#t = t_year - t_year[0]     # option 1 [recommend]: set the first date as zero to avoid float32 precision over tolerance for the interception at the first date\n",
    "t = t_year                 # option 2 also works in Alfredo's data\n",
    "A = np.ones((num_date, 2), dtype=np.float32)\n",
    "A[:, 0] = t\n",
    "vel_est, c = np.dot(np.linalg.pinv(A), dis_sim)\n",
    "dis_est = t * vel_est + c\n",
    "vel_std = np.sqrt(np.sum(np.square(dis_sim - dis_est)) / (np.sum(np.square(t - np.mean(t)))) / (num_date - 2))\n",
    "print('Estimation using scipy.linalg with velocity [m/year]: {:.2f}, velocity std [m/year]: {:.2f}'.format(vel_est, vel_std))\n",
    "\n",
    "## option 2 - scipy.stats.linregress()\n",
    "vfit = stats.linregress(t, dis_sim)\n",
    "vel_est, velstd = vfit[0], vfit[4]\n",
    "print('Estimation using scipy.stats  with velocity [m/year]: {:.2f}, velocity std [m/year]: {:.2f}'.format(vel_est, vel_std))\n",
    "print(vfit)\n",
    "\n",
    "# plot\n",
    "plt.figure()\n",
    "plt.plot(t, dis_sim, '.', ms=8)\n",
    "plt.plot(t, dis_est, '-', lw=4, label='Estimated velocity')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
